#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "util.h"


#define REAL double

REAL* dvector(long nl, long nh)
{
    REAL* v;

    v = (REAL*)malloc((nh - nl + 1 + 1) * sizeof(REAL));
    return v - nl + 1;

}

void free_dvector(REAL *v, long nl, long nh)
{
	free(v + nl - NR_END);
}


REAL** dmatrix(long nrl, long nrh, long ncl, long nch)
{
    REAL** m;
    long i, nrow = nrh - nrl + 1 + 1, ncol = nch - ncl + 1 + 1;

    m = (REAL**)malloc((nrow) * sizeof(REAL*));
    m += 1;
    m -= nrl;

    m[nrl] = (REAL*)malloc((nrow * ncol) * sizeof(REAL));
    m[nrl] += 1;

    m[nrl] -= ncl;

    for (i = nrl + 1; i <= nrh; i++) m[i] = m[i - 1] + ncol;

    return m;
}


void mat2vec(REAL *arr, REAL **mat, int rl, int rr, int cl, int cr)  // 一维数组转二维数组， 一维数组起始索引为0， 二维数组起始终止点[rl, rr]*[cl,cr]
{
	for (int i = rl; i <= rr; i++)
		for (int j = cl; j <= cr; j++)
		{
			arr[(i - rl)*(cr - cl + 1) + (j - cl)] = mat[i][j]; // 一维数组赋值，可以手动计算验证下是否对应，我觉得没问题！！
		}

}

void vec2mat(REAL** mat, REAL *arr, int rl, int rr, int cl, int cr)  // 一维数组转二维数组， 一维数组起始索引为0， 二维数组起始终止点[rl, rr]*[cl,cr]
{
	for (int i = rl; i <= rr; i++)
		for (int j = cl; j <= cr; j++)
		{
			mat[i][j] = arr[(i - rl)*(cr - cl + 1) + (j - cl)]; // 二维数组赋值, 可以手动计算验证下是否对应，我觉得没问题！！
		}

}


void free_dmatrix(REAL** m, long nrl, long nrh, long ncl, long nch)
{
    free(m[nrl] + ncl - 1);
    free(m + nrl - 1);

    return;
}



void mat_add(REAL** a, REAL** b, REAL** c, int xl, int xr, int yl, int yr)
{
    int i, j;

    for (i = xl; i <= xr; i++)
        for (j = yl; j <= yr; j++)
        {
            a[i][j] = b[i][j] + c[i][j];
        }

    return;
}


void mat_add2(REAL** a, REAL** b, REAL** c, REAL** a2, REAL** b2, REAL** c2, int xl, int xr, int yl, int yr)
{
    int i, j;

    for (i = xl; i <= xr; i++)
        for (j = yl; j <= yr; j++)
        {
            a[i][j] = b[i][j] + c[i][j];
            a2[i][j] = b2[i][j] + c2[i][j];
        }

    return;
}


void zero_vector(REAL* a, int xl, int xr)
{
    int i;

    for (i = xl; i <= xr; i++) a[i] = 0.0;

    return;
}

void zero_matrix(REAL** a, int xl, int xr, int yl, int yr)
{
    int i, j;

    for (i = xl; i <= xr; i++)
        for (j = yl; j <= yr; j++)
        {
            a[i][j] = 0.0;
        }

    return;
}

void mat_copy(REAL** a, REAL** b, int xl, int xr, int yl, int yr)
{
    int i, j;

    for (i = xl; i <= xr; i++)
        for (j = yl; j <= yr; j++)
        {
            a[i][j] = b[i][j];
        }
    return;
}

void mat_copy2(REAL** a, REAL** b, REAL** a2, REAL** b2, int xl, int xr, int yl, int yr)
{
    int i, j;

    for (i = xl; i <= xr; i++)
        for (j = yl; j <= yr; j++)
        {
            a[i][j] = b[i][j];
            a2[i][j] = b2[i][j];
        }
    return;
}


void mat_sub(REAL** a, REAL** b, REAL** c,
             int nrl, int nrh, int ncl, int nch)
{
    int i, j;

    for (i = nrl; i <= nrh; i++)
        for (j = ncl; j <= nch; j++)
        {
            a[i][j] = b[i][j] - c[i][j];
        }
    return;
}

void mat_sub2(REAL** a, REAL** b, REAL** c,
              REAL** a2, REAL** b2, REAL** c2,
              int nrl, int nrh, int ncl, int nch)
{
    int i, j;

    for (i = nrl; i <= nrh; i++)
        for (j = ncl; j <= nch; j++)
        {
            a[i][j] = b[i][j] - c[i][j];
            a2[i][j] = b2[i][j] - c2[i][j];
        }

    return;
}

REAL max_matrix(REAL** m, int nxi, int nxt, int nyi, int nyt)
{
    int i, j;
    REAL a = m[1][1];

    for (i = nxi; i <= nxt; i++) {
        for (j = nyi; j <= nyt; j++)
        {
            if (a < m[i][j]) a = m[i][j];
        }
    }
    return a;
}

REAL mat_max(REAL** a, int nrl, int nrh, int ncl, int nch)   
{
    int i, j;
    REAL x = 0.0;

    for (i = nrl; i <= nrh; i++)
        for (j = ncl; j <= nch; j++)
        {
            if (fabs(a[i][j]) > x)
                x = fabs(a[i][j]);
        }
    return x;
}



void print_uvp(REAL** u, REAL** v, REAL** p, int count)
{
    //extern int nx, ny;

    FILE* fu, * fv, * fp;
	
    char buf[20];

    REAL** uu = dmatrix(1, gnx, 1, gny);
    REAL** vv = dmatrix(1, gnx, 1, gny);

    for (int i = 1; i<=gnx; i++)
        for (int j = 1; j <= gnx; j++)
        {
            uu[i][j] = 0.5 * (u[i][j] + u[i - 1][j]);
            vv[i][j] = 0.5 * (v[i][j] + v[i - 1][j-1]);
        }

    sprintf(buf, "data/u%d.m", count);
    fu = fopen(buf, "w");
    sprintf(buf, "data/v%d.m", count);
    fv = fopen(buf, "w");
    sprintf(buf, "data/p%d.m", count);
    fp = fopen(buf, "w");

    print_mat(fu, uu, 1, gnx, 1, gny);
    print_mat(fv, vv, 1, gnx, 1, gny);
    print_mat(fp,  p, 1, gnx, 1, gny);

    free_dmatrix(uu, 1, gnx, 1, gny);
    free_dmatrix(vv, 1, gnx, 1, gny);

    fclose(fu);
    fclose(fv);
    fclose(fp);

    return;
}

void print_mat(FILE* fptr, REAL** a, int nrl, int nrh, int ncl, int nch)
{
    int i, j;

    for (i = nrl; i <= nrh; i++)
    {
        for (j = ncl; j <= nch; j++)
        {
            fprintf(fptr, "  %5.6f", a[i][j]);
        }
        fprintf(fptr, "\n");
    }
    return;
}


void print_bd_posi(REAL** bd_posi, char* bufferposi, int nt, int count)
{
    int k;
    FILE* fbd;
    char buf[20];

    sprintf(buf, "data/bdposi%d.m", count);
    fbd = fopen(buf, "w");

    for (k = 0; k <= nt+1; k++)
    {
         fprintf(fbd, " %f %f \n", bd_posi[k][1], bd_posi[k][2]);
    }
    
    fclose(fbd);
    return;
}


REAL max_vector(REAL* v, int ni, int nt)
{
    int i;
    REAL m = v[nt];

    for (i = ni; i < nt; i++)
    {
        if (m < v[i]) m = v[i];
    }
    return m;
}
